class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 2: EstadÃstica espacial y geoestadÃstica. ### Estimacón de densidad espacial de Kernel (KDE). Gabriel Castro mail gabriel.castro.b@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Septiembre 2023</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ - Problematicas geográficas puntuales. - Estimación de densidad de kernel (KDE). - Ajuste de parametros para la elaboración de un mapa de calor. - Comaparación de bandwidhts. - Kernel con pesos especificos. ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> © Allison Horst ] --- ### Problemáticas geográficas puntuales. #### Análisis de patrones de puntos en R. - Los patrones de puntos son colecciones de puntos geográficos los cuales se asume que han sido generados de forma independiente (Brundson & Comber, 2019). - Algunos de los patrones espaciales más estudiados corresponden a: EpidemiologÃa, riesgos naturales, criminalÃstica, vialidad, permisos de edificación entre otros. - La estimación de densidad de kernel nos permite estudiar patrones espaciales puntuales para identificar zonas con una alta densidad espacial de eventos.  --- ### ¿Que es la estimación de densidad de kernel? - La densidad de kernel (Silverman, 1986) corresponde a una herramienta estadÃstica que calcula la densidad espacial de las entidades considerando los puntos o lÃneas de su vecindad. - El KDE se calcula ponderando las distancias de todos los puntos con una ubicación X,Y en el espacio. Si tenemos más puntos cercanos, la estimación es mayor, lo que indica que existe una alta probabilidad de encontrar un punto en dicha ubicación. - A partir de la distribución de nuestros datos, se crea una "curva suavizada" por medio de un ancho de banda especifico (h) y una función especÃfica del núcleo para la distribución: Gaussiana, cuartica, uniforme, triangular. <center>  <center> <center> ###### (Silverman, 1986) <center> --- ### Tipos de kernel .pull-left[ - Existen distintos tipos de kernel los cuales se diferencian a partir de la función ajustada sobre distribución de nuestros puntos. - Esta función determinará como se va a difundir la influencia de un punto X,Y a lo largo de mi radio de búsqueda. - No existe una función "mejor" o "peor" sino que el ajuste dependerá de la naturaleza de nuestros datos. Hay funciones que se ajustan a zonas más especificas mientras que otras permiten una mayor dispersión de la densidad espacial. - La función "Quartic" o "Biweight" corresponde a la opción más común utilizada para el estudio de patrones espaciales de puntos. (Bailey and Gratrell, 1995; Ratcliffe, 2002; Levine, 2004) ] .pull-right[  <center> ###### (Brian Amberg, 2008) <center> ] --- ### Selección de parámetros para la elaboración de nuestro "mapa de calor". #### KDE requiere de ciertos inputs que son entregados por el usuario los cuales poseen grandes implicancias. - **tamaño de celda**: referido como la resolución de mi KDE. Se encuentra mayormente asociado a la visualización de nuestro mapa, no afecta de sobremanera con los resultados de densidad. Sin embargo, valores muy bajos demandaran mayores recursos en el procesamiento (Chainey, 2013). - **h**: corresponde a nuestro radio de búsqueda (**distancia en unidades de mapa**) que considera una cantidad n de puntos para el "suavizado". Influye de forma directa en la visualización y calculo de la densidad de mi mapa final. - **K / función de ajuste**: corresponde a la función aplicada sobre mi distribución. *Quartic*, *Epanechnikov*, *triweight* etc. --- ### El dilema de nuestro bandwidht. - **h bajos**: Solo los puntos más cercanos a la posición de un punto X,Y recibirán algún peso. Esto dará como resultado una estimación más "ondulada" o "aguda". (Brundson & Comber, 2015) - **h altos**: Esto generara un mayor radio de búsqueda de puntos lo que permitirá que puntos lejanos contribuyan a la representación de densidad espacial. Esto da como resultado estimaciones más chatas y suavizadas. (Brundson & Comber, 2015) <center> <img style="float: width="300" height="360";" src="grafico_bandwith_2.JPG"> <center> --- ### ¿Como seleccionamos nuestro h optimo? #### Existen diferentes métodos: - Selección manual experimental: Probar con diferentes **h** hasta obtener un resultado "visual" optimo (Bailey and Gratrell, 1995 ; Eck et al, 2005) (Conduce a elegir el mapa que se vea mejor) - Selección automática con Rstudio librerÃa MASS: Estima la distancia de búsqueda a partir de un vector de datos (coordenadas). Buen estimador de h para asignar un valor no arbitrario. <b> Desventaja: </b> conduce a valores de **h** sumamente altos llegando a un suavizado muy general. <b> Util </b> para patrones espaciales de puntos muy grandes. - Metodo Brundson & Comber: Estima la distancia de búsqueda a partir de un vector de datos (coordenadas) considerando la desviación estándar y el número total de puntos. <center> <img style="float: width="400" height="250" " src="formula_h_2023.jpg"> <center> --- ### Comencemos a trabajar. #### materiales para la clase. - shapefile accidentes transito Talca. - Shapefile manzanas censales Talca. ```r # Librerias necesarias library(tidyverse) library(sf) library(terra) library(SpatialKDE) library(tidyterra) ``` --- ### Llamar nuestra tabla de sismos y transformarla en un objeto sf. ```r accidentes <- read_sf("accidentes_2020.shp") manzanas <- read_sf("manzanas_talca.shp") ggplot() + geom_sf(data = manzanas) + geom_sf(data = accidentes, colour = "red", size = 1) ```
--- <img src="data:image/png;base64,#DIPGEOPR_02_2_files/figure-html/unnamed-chunk-4-1.png" width="100%" style="display: block; margin: auto;" /> --- ### Definición Bandwidht y creación grilla kernel. ```r # seleccion bandwidht MASS ---------- xcoors <- st_coordinates(accidentes_2021) %>% as.data.frame() h_est <- (bandwidth.nrd(xcoors$Y) + bandwidth.nrd(xcoors$X)) / 2 print(h_est) ``` ``` ## [1] 1419.978 ``` ```r # seleccion Brundson & Comber: ---------- h_x <- sd(xcoors$X)*(2/3*length(xcoors))^1/6 h_y <- sd(xcoors$Y)*(2/3*length(xcoors))^1/6 h_est <- mean(c(h_x, h_y)) print(h_est) ``` ``` ## [1] 327.9913 ``` - Métodos para estimar el h dan resultados muy opuestos a pesar de trabajar con los mismos datos. - La mayorÃa de las funciones presentes en R son para estimar h consideran kernels univariados para una distribución Gaussiana. --- ### Definición Bandwidht y creación grilla kernel. ```r # crear una grilla raster ------- grid_raster <- create_raster(manzanas, cell_size = 100) kde <- kde(accidentes_2021, band_width = 350, kernel = "quartic", grid = grid_raster ) ``` - **Mientras mas pequeño sea el tamaño de la celda, mas tiempo de prosesamiento tendremos.** - Se recomienda mantener las proporciones entre el cell-size de nuestra grilla y el bandwidht. --- ### Visualización del resultado - Debemos transformar nuestro objeto KDE a un spatraster para su visualización con ggplot. - Una vez transformado, reclasificamos como NA todos los valores con densidades cercanas a 0. ```r kde_sprast <- kde %>% rast() kde[kde < 0.000001] <- NA ``` --- ```r p <- ggplot() + geom_sf(data = manzanas, colour = "grey") + geom_spatraster(data = kde, aes(fill = layer), alpha = 0.5, na.rm = T) + scale_fill_whitebox_c("bl_yl_rd") + # Paleta de colores con conversión de valores theme_bw() + labs(fill = "spatial density") + ggtitle("Mapa de calor accidentes de transito Talca 2021") + annotation_north_arrow( location = "tl", # flecha norte which_north = "true", # flecha norte height = unit(1, "cm"), # flecha norte width = unit(1, "cm"), # flecha norte. style = north_arrow_fancy_orienteering( text_col = "black", # flecha norte. line_col = "black", fill = "white" ) ) p ``` --- <img src="data:image/png;base64,#DIPGEOPR_02_2_files/figure-html/unnamed-chunk-11-1.png" width="100%" style="display: block; margin: auto;" /> --- ### Comparemos distintos tipos de kernel. ```r ## ---------------- GRILLA RASTER ---------------------------------------- kde_raster_quartic <- kde(accidentes_2021, band_width = 350, kernel = "quartic", grid = grid_raster) kde_1 <- kde_raster_quartic %>% rast() kde_raster_triweight <- kde(accidentes_2021, band_width = 350, kernel = "triweight", grid = grid_raster) kde_2 <- kde_raster_triweight %>% rast() kde_raster_epanechnikov <- kde(accidentes_2021, band_width = 350, kernel = "epanechnikov", grid = grid_raster) kde_3 <- kde_raster_epanechnikov %>% rast() kde_raster_triangular <- kde(accidentes_2021, band_width = 350, kernel = "triangular", grid = grid_raster, decay = 1) kde_4 <- kde_raster_triangular %>% rast() ``` --- <center> <img src="data:image/png;base64,#Kernels_accidentes.jpeg" alt="drawing" width="1000" height="600"/> <center> --- ### Cambiemos los parametros iniciales - h (bandwidht) - cell-size (tamaño del pixel) ```r grid_raster_H2 <- create_raster(manzanas, cell_size = 30) ## ------------------- CREACIÓN KERNEL ------------------------------------- kde_raster_quartic_h2 <- kde(accidentes_2021, band_width = 200, kernel = "quartic", grid = grid_raster_H2) #--------------------------------------------------------------------------- kde_raster_triweight_h2 <- kde(accidentes_2021, band_width = 200, kernel = "triweight", grid = grid_raster_H2) #--------------------------------------------------------------------------- kde_raster_epanechnikov_h2 <- kde(accidentes_2021, band_width = 200, kernel = "epanechnikov", grid = grid_raster_H2) #--------------------------------------------------------------------------- kde_raster_triangular_h2 <- kde(accidentes_2021, band_width = 200, decay = 1, kernel = "triangular", grid = grid_raster_H2) #---------------------------------------------------------------------- ``` --- <center> <img src="data:image/png;base64,#mapas_h200.jpeg" alt="drawing" width="1000" height="600"/> <center> --- - La la densidad espacial de accidentes disminuyo debido a la disminucion de nuestro h. - Podemos observar las intersecciones con mas accidentes en la ciudad. - Para esta problematica puntual donde es neceario encontrar las zonas mas riesgosas, aplicar un kernel con un suavizado bajo que resalte loa valores centrales es bastante util siendo el **quartic** y **triweight** las opciones mas adecuadas. - Para la escala espacial en la que nos encontramos, la adecuada seleccion de **h** es mas importante respecto a la funcion de kernel que apliquemos. <center> <img src="data:image/png;base64,#3_mapas.PNG" alt="drawing" width="1000" height="350"/> <center> --- ### Densidad de kernel con pesos especificos. - Dentro de la funcion **kde** podemos encontrar el argumento **weights** el cual nos permite generar pesos especificos para cada punto de nuestra área de estudio - Esto resulta de utilidad cuando tenemos problematicas puntuales que poseen algun valor importante a la hora de estudiar su densidad espcial: Magnitud de sismos, personas afectadas por un accidente entre otros. - Para esto el argumento **weights** debe ser de tipo numerico con la misma cantidad de valores que de observaciones. - Columna **Graves** de mi objeto de accidentes para obtener la cantidad de lesionados graves. De esta forma estaremos otorgando un peso mayor dentro de la densidad espacial. --- ### Peso especifico por cantidad de lesionados graves por accidente. ```r # crear una grilla raster ------- acc <- accidentes_2021$Graves acc_graves <- acc print(acc_graves) ``` ``` ## [1] 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [75] 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ## [149] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ## [186] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [223] 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 ## [260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [297] 0 0 0 0 1 1 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ## [334] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [371] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [408] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ## [445] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ## [482] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [519] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [556] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ## [593] 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [630] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [667] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ## [704] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ## [741] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [778] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [815] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [852] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [889] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 ## [926] 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 ## [963] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ## [1000] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ``` --- ```r # crear una grilla raster ------- grid_raster <- create_raster(manzanas, cell_size = 30) # Calculo de KDE con pesos especificos --------------- kde_raster_quartic_weig <- kde(accidentes_2021, scaled = F, band_width = 300, weights = acc_graves, kernel = "quartic", grid = grid_raster ) plot(kde_raster_quartic_weig) ``` --- <center> <img src="data:image/png;base64,#KDE_GRAVES.JPEG" alt="drawing" width="900" height="500"/> <center> - Al tener valores iguales a 0 dentro de los pesos especificos, nuestro KDE nos muestra zonas muy puntuales asociados a la presencia de un valor >0. - Para obetener una densidad mas amplia dentro de la ciudad de Curico, debemos aumentar el tamaño de nuestro **h**. --- ### Cambio de h en nuestro KDE de pesos especificos. - h = 500m <center> <img src="data:image/png;base64,#KDE_ACC_H500.JPEG" alt="drawing" width="900" height="500"/> <center> --- ### BibliografÃa. Chris Brunsdon & Lex Comber, 2015: An introduction to R for Spatial analysis and mapping. Chris Brunsdon & Lex Comber, 2019: An introduction to R for Spatial analysis and mapping. Second edition. Chainey, Spencer, 2013: Examing the influence of cell size and bandwidth size on kernel density estimation crime hotspot for predicting spatial patterns of crime Krisp. J, Peters.S, & Murphy. C, 2009: Visual bandwidth selection for kernel density maps. Xie.Z, Yan.J, 2008: Kernel Density Estimation of traffic accidents in a network space. Density, kernels and occupancy, Source: https://www.spatialanalysisonline.com/HTML/density__kernels_and_occupancy.htm --- class: inverse middle 